A detailed investigation of the properties of a Vlasov-Maxwell 
equilibrium for the force-free Harris sheet 

T.Neukirchfl F. WilsonQ and M. G. Harrisorfl 

School of Mathematics and Statistics, University of St. Andrews, 
St. Andrews KYI 6 9SS, United Kingdom 

Abstract 

A detailed discussion is presented of the Vlasov-Maxwell equilibrium for the force-free Harris 
sheet recently found by Harrison and Neukirch (Phys. Rev. Lett. 102, 135003, 2009). The deriva- 
tion of the distribution function and a discussion of its general properties and their dependence 
on the distribution function parameters will be given. In particular, the distribution function can 
be single-peaked or multi-peaked in two of the velocity components, with possible implications for 
stability. The dependence of the shape of the distribution function on the values of its parameters 
will be investigated and the relation to macroscopic quantities such as the current sheet thickness 
will be discussed. 
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I. INTRODUCTION 



Force-free plasma equilibria with the property 

j x B = — (V x B) x B = 0, (1) 

are of great importance, in particular for space and astrophysical plasmas. Equation ([1]) 
implies that the current density, /ioj = V x B, is parallel to the magnetic field, B, so that it 
can be written as /ioj = ckB. In general, the function a can vary with position, but has to 
be constant along magnetic field lines, since V ■ j = together with V • B = implies that 

B ■ Va = 0. (2) 

Obviously Eq. (j2J) is also satisfied if a = constant. This case is usually referred to as the 
linear force-free case, because the equation determining the magnetic field is linear in this 
case. Magnetic fields for which a varies from field line to field line are called nonlinear 
force-free fields. 

Whereas many force- free equilibria can be found using magnetohydrodynamics (MHD), 
this is not the case when Vlasov- Maxwell (VM) theory is used. Collisionless force-free 
equilibria have only been found for the special case where the magnetic field depends only 
on one spatial Cartesian coordinate (in this paper taken to be z). This case is trivial 
in MHD, but finding the appropriate distribution functions for given magnetic field and 
current density profiles for a collisionless equilibrium is a highly nontrivial task. The reason 
for this difficulty is that one has to try and solve the VM problem in the opposite direction 
than it is usually treated, which is to specify the distribution functions (DFs) and then to 
calculate the magnetic field by solving Ampere's law. 

This difficulty is reflected by the fact that only a very small number of exact force-free VM 
equilibrium DFs are known and all known solutions were of the linear force-free type^ 1 ^ 
until the first nonlinear force-free VM equilibrium DF was presented in a recent Letter.- The 
DFs found in Ref. are for the force-free Harris sheet, with a magnetic shear field ensuring 
force balance instead of a plasma pressure gradient as in the original Harris sheet .— 

For reasons of space no detailed discussion of a) the derivation of the DFs and b) their 
properties has been given in Ref. [fa. In the present publication, we aim to give a full 
discussion of the method used to derive the DFs in Sect. [TT1 and of its properties in Sect. 
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IIHI Of particular interest is the possibility that the DFs can have multiple maxima in two 
of the velocity components (in the coordinate system used in this paper the v x - and v y - 
components), which may have stability implications. Therefore, a detailed investigation of 
the connection between the shape of the distribution function and the parameter values was 
carried out. A summary and conclusion will be presented in Sect. IIVI 



II. CALCULATION OF THE EQUILIBRIUM DISTRIBUTION FUNCTION 



A. Basics 



We use Cartesian coordinates x, y, z complemented by the corresponding velocities v x , 
v y , v z for the DFs. We assume spatial invariance in x and y, i.e. all quantities depend only 
upon z. We also assume time-independence. 

For the problems considered in the present paper the magnetic field has only two non- 
vanishing components, B x and B y , which, using an appropriate gauge, can be written in 
terms of a vector potential A = (A x , A y , 0) in the form 

B. = (3, 

B, = £. (4) 
The electric field is given by the negative gradient of an electric potential <fi such that 

The magnetic and electric fields thus automatically satisfy the homogeneous steady-state 
Maxwell equations V • B = and V x E = 0. 

Due to time independence and spatial symmetries we have three obvious constants of 
motion for particles of species s with charge q s and mass m s moving in these fields, namely 
the particle energy, H s , 

H s = 2 m «( u a + V l + V z) + 1s(t>, (6) 

the canonical momentum in the x-direction, p xs , 

p xs = m s v x + q s A x , (7) 
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and the canonical momentum in the ^-direction, p ys , 

P ys = m s v y + q s A y . (8) 

Solutions of the steady state Vlasov equation 

v.|+*(E + vxB).|=0. (9) 
or m s av 

are given by all positive functions f s depending only on the constants of motion, 

fa = fs(H si p xs ,p ys ), (10) 



and satisfying the appropriate conditions for existence of the velocity moments. If the same 
combination of values for the constants of motion allows particle trajectories in several 
distinct regions of phase space then it is in principle possible to assign different values to 
f s in each region^, but this possibility will not be considered in the present paper (for an 
example of 2D rotationally symmetric VM equilibria see e.g. Ref. 0). 

Using the assumption of quasineutrality to determine the electric potential 0, one can 
show 8 ' 10 that the VM equilibrium problem reduces to solving Ampere's law in the form 



d 2 A x 
d 2 A y 



-Mo 



-Mo 



dP* 
dA, 



where 



Pzz(A x , Ay) — ^ ] TXl s J 



v 2 J s d 3 v 



(11) 
(12) 

(13) 



is the zz-component of the plasma pressure tensor. 

It is obvious (see e.g. Ref LLOj) that Eqs. (ITT]) and (1121) are equivalent to the equations 
of motion of a particle in a 2D conservative potential, with z taking the role of time, A x 
and A y the coordinates of the particle and \1qP zz being the potential. As in the analogous 
particle problem one can integrate Eqs. (TTTj) and (fT2l) once to get 



so that 



d 

dz 



dA, 



2/j, \ dz 



2fi \ dz 



+ 



1 



2/ic 



dA l 
dz 



' J "I" Pzz(A X , Ay) 







(14) 



2^o \ dz 



dA, 



Pzz(A X , Ay) — Pt 



total 



constant, 



(15) 



i.e. the total pressure (magnetic plus plasma pressure) is constant for this class of VM 
equilibria. The total pressure corresponds to the total energy in the particle problem. 

Knowledge of the shape of P ZZ (A X , A y ) can provide insight into the nature of the solutions 
of Eqs. pip and (JT2J) in the same way as knowledge of the potential as a function of position 
in the equivalent particle problem can provide information about the nature of the possible 
trajectories of the particle. It is usually straightforward to calculate P zz as a function of A x 
and A y if the equilibrium DFs are specified. It may, however, also be possible to determine 
equilibrium DFs for a given function P ZZ (A X , A y ) using a method suggested by Channell.- 



B. Channell's Method 



To be able to make analytical progress in determining a distribution function from 
Pzz(A x , A y ) a number of assumptions have to be made. The first assumption made is that 
the dependence of the DFs on the Hamiltonian H s is of the form of a negative exponential, 
i.c 



fs (,H S , P XS , Pyg ) 



nos 



exp(-f3 s H s )g s (Pxs,P ys ) 



(16) 



with (3 S = (ksT s ) 1 , Vth, s = (As?n s ) ^ 2 and g s an unknown function of the canonical mo- 
menta. Using this DF P zz becomes 



Pzz = ~r eM-Psq s (p)N s (A x , A y ) 



with 



N S (A X , Ay) 



n 0s 




exp 



Ps m s / 2 i 2\ 



g S {Pxs,Pys) d,V x dVy 



(17) 



(18) 



The charge density, a, can be calculated by taking the negative derivative of P zz with respect 
to the electric potential-^ as 



a(A x , Ay, (p) = ^2q s exp(-(3 s q s (j))N s (A x , A y 



(19) 



Assuming a two-species plasma with both species having the same charge e with opposite 
sign (e.g. electrons and protons) and quasi-neutrality by letting a = 0, one can determine 
the quasi-neutral electric potential to be 

1 



e(p e + Pi) 



\n(Ni/N e ). 



(20) 



Channell's^ final assumption is strict neutrality, i.e. that Ni(A x ,A y ) = N e (A x , A y ) = 
N(A x ,A y ) for all possible values of A x , A y , implying that <p qn = 0. This will impose 
additional conditions on the parameters of the DFs which have to be satisfied, but this is in 
principle not a problem. 

The neutral P zz is then given by 

PUA x ,A y ) = ^AN(A XJ A y ). (21) 

Using the canonical momenta instead of the velocity components as integration variables 
and using Eq. (l2Tj) . Eq. (TT8l) becomes 



n 0s f f f A 

" — [\Pxs — Qs^x 

PA 




exp | _ ^~-[(^ ~ <lsA x ) 2 + (p ys - q s A y ) 2 ]^ g s (p xs ,p ys ) dp xs dp ys = 

(] e + (3 . P **( A *' A v^ ( 22 ) 

For P ZZ (A X , A y ) a known function of A x and A y , this is a Fredholm integral equation of the 
first type for the unknown function g s (Pxs,Pys)- The kernel K(p xs ,p ys ; q s A x ,q s A y ) of this 
integral equation 

K(p xs ,p ys ;q s A x ,q s A y ) oc exp j _ 7^-[G°zs - q s A x ) 2 + (p ys - q s A y ) 2 ]^ (23) 

depends only upon the difference of its arguments and the standard method for solving such 
integral equations is using Fourier transforms, as also suggested by Channell.- 

It must, however, be pointed out that to be able to determine g s by Fourier transforms 
two conditions need to be satisfied: a) the Fourier transform of P zz (A x ,A y ) must exist and 
b) the inverse Fourier transform to obtain g s must exist. Especially condition b) can prove 
difficult to meet as the inverse Fourier transform involves a factor with the inverse of the 
Gaussian in the convolution integral, i.e. an exponential function with a positive quadratic 
argument. Channel!^ treats several examples for which the Fourier transform method does 
not work using other methods. For the force-free Harris sheet case discussed in the present 
paper, we will also use a more direct method to solve Eq. ( l22l) because Fourier transforms 
are only of limited applicability to our case and because the other method turns out to be 
more instructive. 
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C. Harris sheet and force-free Harris sheet 



The Harris sheet^ is a well-known one-dimensional VM equilibrium. It is widely used 
in theoretical plasma physics, for example for reconnection studies, because it is a typical 
neutral current sheet and is mathematically well-behaved. The magnetic field is given by 

^Harris = £ (tanh(z/lr),0,0), (24) 

the current density by 

^Harris = B O /L(0, 1/ COSh 2 (z/L), 0), (25) 

and the vector potential (in a convenient gauge) by 

A HarHs = B L(0,-ln[cosh(z/L)},0). (26) 

Force balance is maintained by a pressure gradient with P zz (z) given by 

Pq 

Pzz, Harris . o / JT\ Bb,zzi (^^) 

cosh \Zjl) 

with P 0tZZ = 5g/(2/i ) and Pb, zz a constant background pressure. The distribution function 
used by Harris 6 is given by 

fs,Harris = . ^ ~ exp[-(3 s (H s - U ys p ys )], (28) 

{V2irv thiS y 

which is a Maxwellian DF in all velocity directions, but with a constant average bulk flow 
velocity of u ys in the y-direction. Other distribution functions giving rise to the same 



magnetic field and pressure profiles have also been found (see e.g. Ref. QjJ). By using either 



the distribution function ( |28l) directly or Eqs. (1261) and (1271) . one can show that 

zz, Harris 

{A x , A y ) = P , 22 ex V [2A y /(B L)] + P h . zz . (29) 

Note that to get a constant background pressure from the distribution function an extra 
term proportional to exp(— (3 S H S ) has to be added to the right-hand side of Eq. ff28l . Using 
that tanh 2 x — 1 — 1/ cosh 2 x the equilibrium condition (flol) is satisfied with 

B 2 

Ptotal,Harris = 7, ^ Pb,zz- (30) 

Z/i 

The force- free Harris sheet has the same B x as the Harris sheet, but is kept in force balance 
by magnetic pressure due to a B y component, with P zz being constant. The magnetic field 
is then given by 

RffHarris = B G {l&ak{z / L) , 1/ cosh(2:/L), 0). (31) 




FIG. 1: The magnetic field, current density and pressure profiles as functions of z/L for the Harris 
sheet (left panel) and the force- free Harris sheet (right panel). 



The current density is 

HojffHarris = B / L(t&nh(z / L) / cosh(z/L) , 1/ cosh 2 (2/L),0), (32) 
with fioj ff Harris = aB ffHarris where 

= 7 w 7TV - (33) 

L cosh(2/L) 

The vector potential, again in a convenient gauge, is given by 

AffHarris = B L{2 arctan(exp(;z/L) ) , — m(cosh(;z/L) ) , 0) . (34) 

At this point, no form for P zz as a function of A x and A y and no DF are known for this 
equilibrium magnetic field. The derivation of both will be discussed in the next section. 
Plots of the magnetic field components, current density and pressure as functions of z/L are 
shown in Fig. [H 



D. Derivation of the distribution function 



To be able to apply Channell's method to find a DF for the force free Harris sheet, we 
first need to find an appropriate function P ZZ (A X , A y ) for these cases. It can be shown^ that 
to find a P zz that allows a force-free solution is equivalent to finding a potential for which 
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at least one of its equipotential lines is identical with a particle trajectory. The simplest 
examples for this are attractive central potentials whose contours are circles and which also 
allow circular orbits. The corresponding ID VM equilibria are linear force-free magnetic 
fields.i^ 

It is obvious, however, that the P zz for the force-free Harris sheet has to be more complex 



than a central potential. The approach chosen in Ref. 



was to let 



P zz (A x ,Ay) = P 1 (A X ) + P 2 (A y ) 



(35) 



In this case the Eqs. fllip and (fl2l) decouple and can be integrated separately. Thus one can 
see immediately that P 2 (A/) is identical to P zz ,Harris given by Eq. (129]) . The unknown func- 
tion Pi(A x ) can be determined from inverting A x jfHarris(z) using Eq. (134]) and substituting 

Z (A x jf Harris) hltO 



Pi(z) = Pij,- D " 



Using the trigonometric identity 



sin(2a;) 



2/i cosh (z/L) 



2 tanx 



1 + tan x 



2 „i 



one can see that 



sin 



A 



x,ff Harris 

B L 



cosh(z/L) 



so that, dropping the subscript ff Harris, 



P 1 (A x ) = P b>1 -^-sin 2 



A x 



(36) 



(37) 



(38) 



(39) 



2/io \BoL 

Using sin 2 x = [1 — cos(2x)]/2 and putting together Pi(A x ) and P 2 (A y ), we arrive at the 
form of P ZZ (A X , A y ) given in Ref. 



B 2 

Pzz(A X , Ay) = — ~ 



2A n 

- cos , 
2 \B L 



exp 



2Ay_ 

B n L 



Ph 



(40) 



where P = P^i + P 0)ZZ — Bq/{A^q). By construction, Ampere's law (fill and (fl2l) generated 
from this P zz has the vector potential (134|) as a solution, and this solution coincides with a 
contour of P ZZ (A X , A y ). In Fig. [2] we show a surface plot of P ZZ (A X , A y ) for the force-free 
case with the vector potential for the force-free Harris sheet shown as a trajectory at the 
top of the plot. 
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FIG. 2: Surface plot over the A x -A y -plane of the pressure function P ZZ (A X , A y ) for the force-free 
Harris sheet. The vector potential of the force-free Harris sheet traces out a trajectory in the 
A x -Ay-pla,ne, which is shown at the top of the plot. This trajectory coincides with a contour of 
Pzz(A x , A y ), as the general condition for force- free VM equilibria demands*^ 

Having found a P ZZ (A X , A y ), we can use Channell's method 2 to find the corresponding DF. 
As the relation between the unknown function g s {Pxs,Pys) and P zz is linear, it is immediately 
clear that g s must also have the form of a sum, 



g S (Pxs,Pys) = 9sl{Pxs) + 9s2(Pys), 



(41) 



with 



PePi 
Pe + Pi 

Pe + P 



Pi(A x ) 



Ps 



n 0s / exp 



Pa 



2m, 



[Pxs ~ QsA X/ 



9si{Pxs) dp xs 



(42) 



-oo 
oo 



Ps 



2rrm q 



n 0s / exp 



(Pys ~ QsAy) 2 



2m 



QsziPys) dp ys . (43) 



For the time being we can ignore any constant parts of Pi and because the solution for 
a constant P is simply a constant g, which can be added at the end of the calculation due 
to the linearity of the problem. 
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For solving Eq. (j4"2l with P\{A X ) oc cos(2 A x /BqL) one could in principle use Fourier 
transforms, but we shall use a more direct method here. The method is based on the obser- 
vation that, using the trigonometric identity cos[6(s + £)] = cos(fes) cos(fei) — sin(frs) sm(bt), 
one has 

oo 

J exp(— as 2 ) cos[6(s + t)]ds = exp ^— — ^ cos(fet). (44) 

— oo 

Thus, rewriting the integral in Eq. (142]) using v x as integration variable instead of p xs , one 
can see immediately that a g s (p xs ) cos(fl s u xs p xs ) leads to a P\(A X ) oc cos(/3 s Ma; S g s A x ). The 
constant u xs has the dimensions of a velocity so that the argument of the cosine function is 
dimensionless. 

The solution to Eq. fj4*3|) is already known, because this part of the pressure gives rise 
to the y-component of the current density and thus to the Harris sheet B x . Therefore, we 
must have g S 2(p ys ) oc exp(p s u ys pys) (note that the case of a simple exponential P2(A y ) is 
also a special case of one the examples in Channell's paper 2 ). This means that the part of 
the DF depending explicitly on p ys is identical with the p^-dependent part of the original 
Harris sheet DF ([281) . 

The full distribution function therefore has the general form 

f s = r ^ s — — exp(-/3 s F s ) [a s cos {f3 s u xs p xs ) + exp {f3 s u ys p ys ) + b s ] , (45) 
W2-Kv tKs ) 6 



with a s , b s , u xs and u ys being constant parameters of the DF in addition to n 0s and /3 S . We 
remark that we assume that b s > \a s \ > at this point to ensure that f s remains positive. 
The parameters of the DF will have to satisfy a number of consistency relations due to the 
assumptions made for applying Channell's method and in order to relate the microscopic 
DF parameters to the macroscopic parameters B and L. 

E. Consistency Relations 

The pressure tensor component P zz we obtain using Eq. fT4o^ is of the general form f[T7|) 
with 



N s (A x ,A y ) = n 0s exp 



f3 s m s u ys 
2 



2 N r ( P s m s (u 2 xs + u: y 

a s exp — cos(p s u xs q s A 2 



2 



+ exp((3 s u ys q s A y ) + b s exp 
11 



Psm s u 2 ys 



(46) 



The fundamental condition for Channell's method to be applicable is N e (A x , A y ) = 
Ni(A x , A y ). This is satisfied if 

n 0e exp I - 1 = n 0i exp I - 1 = n (47) 

( Pem e (ul e + u 2 )\ ( Am 4 (t4 + ^)\ 
a e exp I — = a« exp I — = a (48) 



2 J V 2 
b e exp — -2- = 6, exp — = 6 (49) 



2 

/3eKe| = (50) 
Pe u ye = ftiUyi (51) 

For the case of the original Harris sheet, Eq. (|5T!) is well known^ as the condition for 
a vanishing electric potential. In the Harris sheet the constant average bulk 

velocity of species s in the ^-direction and condition (15T1) is basically specifying a particular 
frame of reference. In the case of the force-free Harris sheet, the average bulk velocity for 
both the x- and the y- velocity components varies with z and one thus needs more conditions, 
but in principle one can still interpret Eqs. (|47|) to (15T|) as conditions for a particular frame 
of reference in which the electric potential vanishes. 

Using Eqs. ( 1471) to ( 15T1) the general expression for P ZZ (A X , A y ) for the force-free Harris 
sheet equilibrium becomes 

P ZZ (A X , A y ) = — ''J' f f % n [a cos(ep e u xe A x ) + exp(-ef3 e u ye A y ) + b] , (52) 

where, for simplicity, we have used the electron parameters only at the moment. An ex- 
pression which is symmetrical in the electron and ion parameters will be derived in Sect. 

ME 



III. PROPERTIES OF THE EQUILIBRIUM DISTRIBUTION FUNCTION 

A. Relation between microscopic and macroscopic parameters 

Although we have now derived the DF for the force- free Harris sheet, we have not yet 
related the set of microscopic parameters of the DF, namely /3 S , u xs , u ys , a s and b s , to the 
macroscopic parameters of the equilibria, which are B and L. The easiest way to find this 
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connection is to compare Eq. (|52|) with Eq. (|40l) . This leads to 

Bl . ft+ft 



2 W Aft (53) 

i = o, (54) 

A = %*n A (55) 

2 

e/3 e |w xe | = e(3i\u xi \, (56) 



|5 |L 

e(3 e Uy e = ePiUyi, (57) 



2 



5 L 

where we have assumed that L is positive, but allow for Bq to be negative. 

To make the connection with the original Harris sheet results we use Eqs. ( )57|) and (|53|) 



to derive an expression for L in the form (see also Ref. [l2j, Chapter 6) 



L=[ ..*^ M . S \ (58) 



fi e 2 /3 e f3in (uyi - u 



ye I 



which is symmetric in electron and ion parameters. Using Eq. (l57j) . expressions for L using 
only electron or only ion parameters can be derived from Eq. ( 1581) . 

The relation of the other macroscopic parameters to the microscopic parameters are more 
obvious. Equation (1531) directly relates B , the magnetic field strength in the limit z — > oo, 
with (3 e , fa and n . In the original Harris sheet case, n is the maximum value of the in- 
dependent part of the particle density at z = 0, and Eq. (|53|) simply states that the magnetic 
pressure for z —>■ oo has to be equal to the plasma pressure at z = due to force balance. As 
we will see later, in the force-free Harris sheet case the meaning of no changes, but because 
we have effectively separated the total force-balance into two conditions for B x and B y , the 
same condition as for the original Harris sheet still applies for the force-free Harris sheet as 
well. 

Equation fl54|) directly shows that for the force-free Harris sheet we have a = 1/2. 

Equation (1551) relates the constant background pressure P& to the microscopic parameter 
b, which is representing the magnitude of the part of the DF which depends only on H s . Ob- 
viously, b is simply the ratio of the background pressure P& to the pressure (/3 e + Pi) n o/ (A:A)- 
Equation (|56|) . together with Eq. (157|) . allows us to relate u ys to u xs by writing 

\Uys\ = \Uxs\- (59) 



13 



An expression for N(A X , A y ) for the force- free Harris sheet equilibrium which is symmet- 
rical in electron and ion parameters is given by 



N(A x ,A y 



with 



no 



A 



a cos 



2A q 



A 

2(& + A) 



+ exp 



2 Ay 

A 



(60) 



(61) 



An expression for P zz which is symmetrical in ion and electron parameters is obtained by 
using (I6"U|) in Eq. (T2TT) . Using the vector potential for the force-free Harris sheet, ( 1511) . we 
obtain for the particle density as a function of z, expressed using microscopic parameters, 



N(z) = n e (z) = rii{z) = n 



5" 



(62) 



the pressure P zz is obtained by multiplying N(z) by (/3 e + j3i)//3 e [3i. 

The mean bulk flow velocities of each species in the x- and y-directions as functions of z, 
namely 

u ys sinh(z/L) 



< v T , > 



< Vy S > 



;i + 6) cosh 2 (^/L) 



u 



[\ + b) cosh 2 (^/L) 



(63) 
(64) 



which gives a current density of the form 



en {u yi - u yej 



sinh(z/L) 
cosh 2 (z/L) 
1 



J y = en {u yi - Uy 



(65) 
(66) 



ye cosh 2 {z/LY 

The force-free parameter a(z) can be directly determined by using (|58|) in (|33|) resulting 



in 



a(z) 



Hoe 2 (3 e f3jn (uyi - u ye ) 
2(J3 e + /3i) 



2\ 1/2 



cosh 



fj j0 e 2 p e p i n (uyi - u ye ) s 
2(& + A) 



2\ 1/2 



(67) 



One can easily show that this is consistent with the expression for a(z) obtained from the 
current density fl65l) and (166]) and the magnetic field for the force- free Harris sheet fl3Tl) . 
when taking into account ( 15*31) . 
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B. The number of maxima of the DF in v x and v y 

One of the interesting features of the force-free Harris sheet DF (|45j) is that it can have 
multiple maxima in both the v x - and the indirections. We shall discuss the indirection 
first as it is simpler to understand. Looking at the structure of the DF in the v y direction 
one can immediately see that it consists of the Harris sheet DF part, which is a Maxwellian 
distribution function drifting with a constant velocity u ys in the f y -direction, and a part 
which, if regarded purely as function of v y , is Maxwellian at rest. It is intuitively clear that 
one should get two maxima in v y if the drift velocity u ys increases, because the drifting 
Maxwellian moves towards the tail of the Maxwellian at rest. As we show in appendix [A] 
it is relatively straightforward to work out that a necessary condition for having more than 
one maximum in the f^-direction is 



i.e. the constant drift velocity has to be larger than twice the thermal velocity. There is, 
however, a second condition on the parameter b s that also needs to be satisfied for the DF 
to have more than one maximum in v y . We derive and state the exact condition in appendix 
|A~1 but its physical meaning is very easy to understand. If h s exceeds a certain limiting value, 
the part of the DF which does have vanishing average velocity in the ^-direction dominates 
over the other part of the DF, so that a second maximum cannot develop even if (1681) is 
satisfied. Usually, this condition on b s will not be very restrictive, though, as the upper limit 
for b s grows exponentially with Uy S /vf hs (see appendix [A}. We show examples of DFs as 
functions of v y for the different cases in Figs. [3] - [51 For these figures the values of b s have 
been chosen to be close to the critical value discussed in appendix [A] for illustrative purposes. 
The values for b s are of the order 4.5 • 10 3 for the examples shown, which corroborates the 
point made above regarding the exponential growth of the limiting value. 

We now turn to the dependence of the DF on v x . Due to the cosine-dependence it is clear 
that the possibility of multiple maxima in v x exists. We discuss the details of the calculation 
in appendix [HI From the analysis in appendix [B] we find that the condition for having just 
a single maximum in v x is 



u 



> 2v, 



th,si 



(68) 




(69) 
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FIG. 3: Shape of the DF in the Vy-direction for various values of z/L for a multiple maximum case. 
Here u ys = 2>vth,s: b s = 4.254 • 10 3 and v x = have been used. In the case shown here the DF has 
multiple maxima for small values of \z\, but only a single maximum as \z\ increases. The value of 
b s has been chosen to be smaller than, but close to the critical value calculated in appendix [Al 

This condition on b s can be understood in the same way as the similar condition on b s 
derived for the ^-dependence. If b s is large enough the Maxwellian background plasma it 
represents dominates the part of the DF with the cosine dependence and we only have a 
single maximum of the distribution function. If the condition (I69I) is not satisfied then we 
have multiple maxima in v x , but their existence still depends on the values of z/L and v y . 
Obviously, for small u ys /v t h, s the limiting value on the right hand side of (1691) is 1/2, which 
is consistent with the absolute lower limit on b s mentioned before. Examples of the different 
cases are shown in Figs. [6]to[8j 

A slightly different perspective on the discussion above can be provided if we express the 
ratio Uy S /v t h,s in terms of the current sheet thickness L. Using Eq. (l57l) . we get 

u y s — a r ^ s 



4-£, (70) 



where r 9iS = m s v t h,s/eB is the thermal gyroradius of species s. If all parameters except L 
and u ys are fixed, it is obvious that a decrease in the current sheet thickness will eventually 
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FIG. 4: Shape of the DF in the fy-direction for various values of z/L for the critical case, at 
which the transition between multiple maxima and a single maximum occurs. Here u ys = 3i>th,s> 
b s = 4.427 • 10 3 and v x = have been used. For z = the DF has a point of inflection with 
horizontal slope, but only one maximum. For \z\ > the DF only has a single maximum. The 
value of b s has been chosen to be equal to the critical value calculated in appendix 1X1 

lead to multiple maxima in the DF, first in v x by violating condition fl69l) and then in v y as 
well. This may obviously have implications for possible velocity instabilities of the system, 
e.g. the two-stream or bump-on-tail instabilities, apart from macroscopic instabilities of the 
current sheet, e.g. the collisionless tearing mode. A detailed investigation of the stability 
properties of this inhomogeneous Vlasov- Maxwell equilibrium would be very interesting, but 
is beyond the scope of the present paper and will be left for future work. 



IV. SUMMARY AND CONCLUSIONS 



We have given a detailed presentation of the derivation and the properties of the DF 
for the collisionless force-free Harris sheet found in Ref. 0. In particular, we have shown 
how the microscopic parameters of the DF are related to the macroscopic parameters of the 
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FIG. 5: Shape of the DF in the Vy-direction for various values of z/L for a single maximum case. 
Here u ys = 3vth,s, b s = 4.659 • 10 3 and v x = have been used. The value of b s has been chosen to 
be greater than the critical value calculated in appendix lAl 



z = z = 0.5 
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FIG. 6: Shape of the DF in the t^-direction for various values of z/L for a single maximum case. 
Here u ys = vth,s, b s = 2.85 and v y = have been used. 
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FIG. 7: Shape of the DF in the ^-direction for various values of z/L for a multiple maximum 
case. Here u ys = vth,s, b s = 1.43 and v y = have been used. In the case shown here the DF has 
multiple maxima close to the sheet centre {z = 0), but a single maximum for larger distances from 
the sheet centre. 
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FIG. 8: Shape of the DF in the w z -direction for various values of z/L for a multiple maximum 
case. Here u ys = 2vth,s, b s = 28.66 and v y = have been used. In the case shown here the DF has 
multiple maxima for all values of z. 
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magnetic field. We have also given a detailed derivation of the conditions on the parameters 
of the DF to ensure that it has only a single maximum in v x and in v y . We have shown that 
as the current sheet thickness decreases the condition for multiple maxima will eventually be 
violated and we have suggested that this may lead to velocity space instabilities in addition 
to other macroscopic instabilities for thin current sheets. The stability properties of the VM 
equilibrium are a very interesting topic for further investigations. 
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APPENDIX A: CONDITION FOR TWO MAXIMA IN THE ^-DIRECTION 



From a mathematical point of view it is easier to write the DF as a function of the 
momenta when carrying out this calculation. The p 2S -dependence does not play any role in 
the calculation and can be integrated out. The reduced DF for p xs and p ys then reads 

1 



F s (z,p xs ,p ys ) = exp 



(Pxs A x 

) + {Pys ~ A/) | \ i a s cos (Pxs) + exp(p ys ) + b s ] , 

(Al) 



where F s = 2ir(m s v t h,s) 2 F s /n 0s , with F s = j f s dp 



zsi Uys 



^ys I ^th,si Pxs PsUysPxsi Pys 



PsUysPys, A x = q s (3 s u ys A x = 2A X /(B L) and similarly A y = q s (3 s u ys A y = 2A y /(B L). 
For an extremum of F s in the p ys direction the derivative 



= exp 



1 



(PXS Axj + (p yS Ay) 

1 



X 



exp(p ys ) - (p ys - A y ) [a s cos(p xs ) + exp(p ys ) + b s 

U ys 



(A2) 



must vanish, leading to the condition 

Pys ~ Ay = 



u 2 ys exp(py S 



a s cos(p xs ) + exp(p ys ) + b s 
We remark that the right hand side is well-defined because b s > a s > 0. 



(A3) 
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The left hand side of ( 1A3I) is a linear function of unit slope in p ys , which crosses the 
p ys -axis at p ys = A y . As A y varies between — oo and 0, the left hand side intercepts the 
p ys -a.xis for negative values of p ys . The right hand side of (1A3|) can be rewritten as 

= H ^) ' <A4) 
where A = u ys > and B = a s cos(p xs ) + b s > 0. The function (IA4I) is positive, increases 
monotonically and is bounded between and A. Therefore, a necessary condition for mul- 
tiple maxima of the DF in v y (or p ys ) is that the maximum slope of R{p ys ) must be larger 
than 1. Otherwise the flA3j) can only have a single solution, implying a single maximum 
for the distribution function. It is straightforward to show that R{p ys ) has its maximum 
slope A/4 at p ys = In B. So the necessary condition for multiple maxima is A/A > 1 which 
translates into 

\u ys \ > 2v th>s . (A5) 

However, (1A5I) is not sufficient, because even if it is satisfied, it is still possible that 
p ys — A y intersects with R{p ys ) only once, namely if the value of B is large enough. As 
discussed above the left hand side of (IA3I) can only cross the for p ys < 0, depending 

on the value of A y (and thus z/L). Since the In B is positive it can happen that R{p ys ) takes 
on its maximum slope too far to the right for more than one intersection between the two 
functions to happen. The transition between three intersections to one intersection happens 
at the value of B for which the straight line of slope one through the origin just touches 
the graph of R(p ys ) at the point where it also has unit slope (see Fig. [9]). One can easily 
calculate the value of p ys for which the function R{p ys ) has unit slope as 



p ys , u = \n(2B) - ln(A - 2 - y/A(A-A)). (A6) 
Two remarks are to be made here: 

• Pys,u only has a real value if A > 4, which is consistent with the condition found before 
for R(p ys ) to have slope greater than unity anywhere; 

• For A > 4, the function R{p ys ) has unit slope at two values of p ys , of which one has 
to choose the larger one (see Fig. [9]), as we have done above. 

The limiting value for B can now derived from 

Pys,u R(pys,u) (A7) 
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FIG. 9: Upper left panel: A case for which R(p ys ) has a maximum slope smaller than one (A = 2.0, 
B = 2.0). Upper right panel: A case for which R(p ys ) has a maximum slope equal to one (A = 4.0, 
B = 2.0). Lower left panel: A case for which R(p ys ) has a maximum slope larger than one (A = 6.0, 
B = 2.0). A case for which R{p ys ) has a maximum slope greater than one, but for which B is larger 
than B\ (A = 6.0, B = 33.5, B\ = 30.42). The straight line shown passes through the point of 
maximum slope in all plots apart from the lower right panel. In the lower right panel the straight 
line passing through the origin is shown. 



leading to 



Bi = \[A - 2 - v^-4)] exp ^ 



2A 



A-^A(A-A) ' 



(A8) 



so the sought for condition is 



B<Bl (A9) 
Since B still depends upon p xs we have to replace it by the minimum value it can take on 
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as function of p xs to get a condition which is independent of p xs . 

In summary, the DF has more than one maximum in p ys (and thus in v y ) if the following 
conditions are both satisfied 

\Uy S \ > 2V t h,S, ( Al0 ) 



2 



b s < 7T—( U l S - 2v th,s - \ U Vs\\l U l S - Av lh, 8 ) eX P ' — ~ 



+ 5 e X p^J, (AH) 
where we have made use of fljSP to replace a s 



APPENDIX B: CONDITION FOR MULTIPLE MAXIMA IN THE 
DIRECTION 

The analysis here is very similar to that in appendix [A] Again we use the reduced DF 

(I A II) expressed as a function of the canonical momenta p xs and p ys . Taking the derivative 

of F s with respect to p xs gives 

dF s f 1 r , T x 2 /_ x x 2] 1 

^- = -ex P j-— [(p„. - + ( Pys -A y )\jx 

| a s sin(p xs ) + i (p xa - A^) [a s cos(p a . s ) + exp(p ys ) + b a ] \ . (Bl) 

I u ys J 



Setting this to zero gives the equation 



u 2 ys a s sm{p xs 



p xs ~A X = "° 7-\ + h > ^ 

a s cos(p xs ) + exp(p ys ) + b s 

or, in an abbreviated form 

p xs - A x = R(p xs ), (B3) 

with 

R(P*s) = -^p^L (B4) 
yPxsJ cos(p xs ) + D 1 ; 



where C = u ys > and D = (b s + exp(p ys ))/a s > 1, because b s > a s . The function (JB4|) i 



is 



a bounded periodic function of p xs . Furthermore, A x = 4 arctan(exp(z/L)) varies between 
and 2-7T, so the left hand side of (1B3[) can only cross the p xs -axis between and 2tt. The 
slope of R(p xs ) is given by 

dR_ = Dcos(p xs ) + 1 
(cos(p xs ) + D) 2 ' 
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FIG. 10: Left panel: A case in which R{p xs ) has a maximum slope greater than unity (A = 1.0, 
B = 1.5); Middle panel: The limiting case with maximum slope equal unity [A = 1.0, B = 2.0); 
Right panel: A case in which R(p xs ) has a maximum slope less than unity (A = 1.0, B = 2.5). For 
these plots the straight line of unit slope has been chosen to cross the p^-axis at p x 



7T. 



which shows that R{p xs ) has a positive slope for cos{p xs ) < —1/D, which is always satisfied 
for some p xs in the interval < p xs < 2n. Therefore, a necessary and sufficient condition for 
multiple maxima of the DF in v x is that R(p xs ) has a maximum slope which is larger than 
unity. Examples for the different cases are shown in Fig. [TUJ 
Taking the derivative of (1B5[) we get 

— ~ C zm(p xs ) D2 - 2 - D cos fe*) (B6) 
dp 2 xs xs (cos(p xs ) + D) 3 

A brief calculation shows that R(p xs ) has positive slope only for p xs = mr with n an odd 
integer. The maximum value of the slope is given by C/(D — 1), which leads to the condition 



C < D- 1 



(B7) 



for the DF to have only one maximum. The lowest value D can take (as a function of p ys 
is D = b s /a s so that we finally arrive at the condition 



b s > — exp 




(B8) 

for the DF to have only one maximum in v x , where we have used ( I48p and fl54l) to replace 
a,. 
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